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ABSTRACT 

Dark matter-dominated cluster-scale halos act as an important cosmological probe and provide a key testing 
ground for structure formation theory. Focusing on their mass profiles, we have carried out (gravity-only) simu- 
lations of the concordance ACDM cosmology, covering a mass range of 2 x 10 12 -2 x 1O 15 /2 _1 M and a redshift 
range of z = 0-2, while satisfying the associated requirements of resolution and statistical control. When fitting 
to the Navarro-Frenk- White profile, our concentration-mass (c—M) relation differs in normalization and shape in 
comparison to previous studies that have limited statistics in the upper end of the mass range. We show that the 
flattening of the c-M relation with redshift is naturally expressed if c is viewed as a function of the peak height 
parameter, v. Unlike the c-M relation, the slope of the c-v relation is effectively constant over the redshift 
range z = 0-2, while the amplitude varies by ~ 30% for massive clusters. This relation is, however, not universal: 
Using a simulation suite covering the allowed wCDM parameter space, we show that the c-v relation varies by 
about ± 20% as cosmological parameters are varied. At fixed mass, the c(M) distribution is well-fit by a Gaussian 
with a c /(c) ~ 0.33, independent of the radius at which the concentration is defined, the halo dynamical state, and 
the underlying cosmology. We compare the ACDM predictions with observations of halo concentrations from 
strong lensing, weak lensing, galaxy kinematics, and X-ray data, finding good agreement for massive clusters 
(M V! > > 4 x 10 14 /! _1 M ), but with some disagreements at lower masses. Because of uncertainty in observational 
systematics and modeling of baryonic physics, the significance of these discrepancies remains unclear. 

Subject headings: Cosmology: clusters-profiles — methods: A^-body simulations 



1. INTRODUCTION 

According to the current cosmological model, structure forms 
in the Universe primarily by the amplification of primordial 
fluctuations driven by the gravitational Jeans instability. The 
process of nonlinear structure formation is hierarchical and com- 
plex, the initial perturbations evolving eventually into a 'cosmic 
web' network consisting of voids, filaments, and clumps. The 
clumps, termed halos in cosmological parlance, are dark matter 
dominated localized mass overdensities with their own com- 
plex substructure. Observed baryonic systems such as galax- 
ies and hot gas reside in these halos. Although the dark mat- 
ter within halos cannot be observed directly, its presence can 
be inferred by dynamical arguments, and much more directly, 
through gravitational lensing of background sources. 

The notion of the dark matter dominated halo is one of the 
fundamental building blocks in studies of the formation of in- 
dividual gala xies, galaxy groups, and galaxy clu sters (for an 
overview, see lMo. van den Bosch. & White Il2010l) . The struc- 
ture of halos has been extensively studied using N-body simula- 
tions over a wide range of halo masses. Even though individual 
halos can be, and are, dynamically and morphologically com- 
plex, it was shown by iNavarro. Frenk. & White 1 1 19961 1 19971) 
(NFW) that the spherically averaged density profiles of 're- 
laxed' halos formed in cold dark matter (CDM) simulations can 
be described by a roughly universal functional form - the NFW 
profile - independent of their mass, the spectrum of initial fluc- 
tuations, and cosmological parameters. The NFW profile has 
a fixed shape, albeit with two scale parameters; as applied to 
individual halos it has been remarkably successful and is often 
applied to all halos, regardless of their dynamical state. (When 
applied to stacked or average halos, this profile is somewhat 



less succesful, as discussed later below.) 

The two parameters of the NFW profile are a mass and a scale 
radius. The scale radius, r s , specifies the point where the loga- 
rithmic slope of the profile equals -2 (at small radii, the profile 
~ l/r, while at large radii, it asymptotes to ~ 1 /V 3 ). Instead of 
r s , one often uses the concentration, which is the radial scale set 
by the halo mass divided by r s . In cluster cosmology, the usual 
key observable is the halo mass, rather than the profile per se. 
The cluster mass function (cluster abundance, more generally), 
is a sensitive probe of dark energy, since clusters form very late, 
during the epoch of dark energy dominance. However, measur- 
ing the concentration parameter, the simplest first measurement 
of a profile, can also be very useful. 

First, as shown originally by NFW, the concentration of a 
halo, c, has a strong correlation to its mass, M, therefore mea- 
suring the c-M relation observationally is a direct test of the 
CDM paradigm. In fact, combining cluster c-M predictions 
and measurements, and the m easured gas mass f raction, one 
can aim to constrain VL m and cr% dEttori et al. 11201 ll) . As another 
example, lensing shear peak counts, a proposed weak lensing 
survey cosmol ogical probe, is very s ensitive to the form of the 
c-M relation dKing & Mead 1120111) . Finally, future measure- 
ments of the weak lensing power spectrum will probe small 
enough spatial scales that results will be sensitive to baryonic 
effects on the halo profile, i.e., modifications to the gravity-only 
c-M relation dWhite ll2004tlZhan & Knoxll2004 . We will re- 
turn to these points in more detail below. 

The correlation of halo concentration with mass is based on 
the idea - as first explicated by NFW - that the concentration 
is determined by the mean density of the universe when the 
halo is assembled, with higher concentrations corresponding to 
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higher densities. Thus cluster mass halos, which are assembling 
today, should have a lower concentration than halos of lower 
mass that were built up at an earlier epoch, where the mean 
density was higher. Furthermore, one may expect this trend to 
flatten out (sufficiently) beyond the nonlinear mass scale M*, 
and therefore, since falls rapidly with redshift, flatten out 
over an extended range in mass as redshift increases. Although 
the general arguments are plausible and are broadly consis- 
tent with simulation results, a predictive theory for the mean 
of the c-M relation, and its scatter, does not exist. Several sim- 
ple he uristic models tuned to simulations have been suggested 
(NFW; iBullocket al.ll200U lEke et al. Il200lt IZhao et al Il2009h 
but their predictive status cannot be considered satis factory, es 
pecially at the higher end of halo masses (see, e.g., Gao et al 



20081; lHavashi & White 120081; iMaccio. Dutton, & van den Boschl 
2008t IZhao et al 1120091) . Indeed there is sufficient uncertainty 



even when comparing simulation results from different groups, 
that the general problem is still open. However, as the mass 
resolution in large-volume simulations continues to improve, 
we may expect this situation to be merely temporary. 

On the observational front, cluster (and group) halo profiles 



individually, and in combination (see, e.g., Comerford & 


. Nataraianl 


2007 


; Broadhurst et alj2008; Mandelbaum et al. 12008; 


Okabe et al. 


201C 


;Oeurietal. 12011; Zitrin et al. 201 it Coe et al. 1 2012), 



projected gas densit y and temperature prof i les from X-ray ob- 
servations (see, e.g.. IVikhlinin et al. Il2006t iBuote et al 20071; 



Schmidt & Allen 120071; iGastaldello et al 120071; I Vi khlinin et al. I 
20091; ISun et al. ll2009HEttori et al. 1201 lb, and galaxy kinemat- 
ics dDiaferio, Geller, & Rines 1 120051; iRines & Diaferio I l2006t 
IWoitak&Lokasll2010l and references therein). Results from 
these observations have generally shown qualitative agreement 
with the c-M relation obtained from simulations, although 
there have been difficulties with matching the shape and nor- 
malization. Additionally, there are discrepancies between dif- 
ferent sets of observations, presumably because the underlying 
systematics are not fully understood and modeled. 

The purpose of this paper is to present a set of predictions 
for the NFW mass profile targeted primarily towards massive 
clusters. To do so, however, a fairly large mass range must be 
considered in order to obtain a sufficiently well-defined c-M 
relation. Our simulations cover three orders of magnitude in 
mass (~ 10 12 - ~ lO 15 /2 -1 M ) with very good control of statis- 
tics over the entire range. The high dynamic range and excellent 
statistics enable us to derive a new set of results for the mass 
profile, including profile evolution and probability distribution 
functions (PDFs) for the concentration as a function of mass. 
We compare our results with previous simulations and with a 
set of recent observations of the cluster mass profile. 

The paper is organized as follows. In Section [2] we discuss 
general features of the c-M relation in the simulation context 
focusing on the role of differing definitions and analyses. In 
Section[3] we describe the main features of the simulation runs. 
We present our results for the c-M relation and its redshift 
evolution in Section [4] This is followed (Section [5]l by a pre- 
sentation of results from a suite of wCDM cosmologies in order 
to further study how the concentration depends on cosmology. 
Next, in Section [6] we provide a detailed comparison with re- 
cent observations, noting areas of agreement and disagreement. 
Finally, Section [7] is devoted to a summary of the results and 
further discussion. An Appendix discusses various systematic 
issues that need to be considered when deriving concentrations 



from simulation results. A number of tests are used to illus- 
trate these points and to verify the robustness of the numerical 
procedures carried out in this paper. 

2. HALOS AND CONCENTRATIONS 

Dark matter dominated halos are dynamically complicated 
and rendering them as simplified 'few parameter' objects in- 
volves a fair degree of approximation, opening the possibility 
of biases in the sense that different procedures will inevitably 
yield different results - what these different results may imply 
for observations is yet another question. In this paper we adopt 
a minimal approach to describing halos; we consider the first 
approximate description of a halo to be a spherically averaged 
profile with a single power law and one overall parameter (e.g., 
singular isothermal sphere), and the NFW profile as essentially 
taking the next step with a broken power law and two param- 
eters (the mass and the concentration). In three dimensions, 
halos are known to be triaxial with a major axis r oughly twice 
as lon g as the two minor axes (roughly equal) dJing & Suto I 
120021) . Spherical averaging of this profile yields the NFW bro- 
ken power-law. 

In reality, halos have complicated substructure and complex 
infall regions, all of which may make interpreting the concen- 
tration somewhat nontrivial, as well as introduce projection- 
related biases in observation s (e.g., White. Cohn., & Smit |2010|) . 
Nevertheless, as shown by lEvrard et al. I ( 120081) . cluster-scale 
systems with mass greater than 10 14 /! _1 M© are dominated by 
large, primary halos - satellite halos carrying only ~ 10% of 
the mass - and possess a well-defined and regular virial rela- 
tion. Therefore, it appears reasonable to proceed in the manner 
outlined above. 

The lack of smoothness in the individual radial density pro- 
files of halos - even at high mass resolution - means that the 
simple N FW description will have varying levels of success 
(see, e.gjTormen, Bouchet & Whitell997tlLukic et al. Il2009t 
iReed. Koushiappas. & Gao 11201 ll) on a halo by halo basis. Av- 
erage or stacked profiles are of course much smoother; it turns 
out that such profiles systematically deviate from the NFW form 
and another scale parameter is often introduced to i mprove the 
fit, le ading to the so-called Einasto profile (see, e.g jGao et al. I 
120081) . While this improves the stacked fit primarily at smaller 
radii, it has little ef fect on measurements of the concentration 
for in dividual halos dGao et al. l2008l ; lReed. Koushiappas. & Gao I 

Since our objective is to carry out comparisons primarily 
against observations of individual objects, rather than against 
correlation functions or stacked observations, we do not use the 
Einasto profile. 

An important piece of missing physics in our simulations is 
the lack of non-gravitational baryonic effects. This is a very 
difficult problem to deal with for galaxy and group-scale ob- 
jects, but less so for clusters. In clusters, the dominant form 
of atomic matter is not stars, but hot gas. Gas cooling does 
not have a major effect on the profiles excep t close to the inner 
regions of the cluster, roughly r < 0.1R v / r dKazantzidis et al.l 
l2004t iDuffv etaT~ll2010l: ICui et al. II20TTI) . Beyond this radius 
the gas distributio n is determined by t he self-consistent gravi- 
tational potential. IDuffv et al. I d2010l) have carried out an ex- 
tensive study of possible baryonic effects (cooling, feedback) 
on cluster profiles and concluded that the baryonic effects are 
likely to alter the concentration at most at the 10% level. This is 
roughly the level of systematic control over the current gravity- 
only measurements of halo concentrations, therefore we do not 
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concern ourselves with estimating baryonic effects or trying 
to correct for them, beyond not fitting for the concentration at 
radii, r < 0.1/? V! >. The fact that we have good agreement with 
observations for massive clusters (Section [6]) may be viewed 
as added support to the argument that baryonic effects do not 
influence cluster profiles away from the inner regions. 

A large number of numerical studies have been carried out 
investigating halo profiles and paying close attention to the be- 
havior of the density cusp on the very smallest scales. We are, 
however, concerned not with these scales, but more with scales 
of order ~ (0.1- l)R v t r , since our target halos have relatively 
modest concentrations. Previous numerical simulations have 
found that in the region of interest to us, the concentration is a 
slowly varying function of mass, typically described by power 
laws with index a ~ -0. 1 at z = 0. These simulations have var- 
ied widely in dynamic range, box size, and mass resolution. 
Partly as a result of this, there have been some disagreements 
in the value of the slope and the normalization of the c—M 
relation, and also some lack of clarity regarding the reasons un- 
derlying the differences. 

Among the more rece nt studies are those in volving the Mil- 
lenium simulation (MS) dSpringel et al. 1120051) with 2160 3 par- 
ticles an d a box of side 500 /r'Mpc assuming a WMAP1 cos- 
molo gy (iNeto et al. Il2007t iGao et al. Il2008t lHavashi & White I 
2008). Halo profiles were investigated over a mass range of 
1O 12 -1O 15 /i _1 M and it was found that a ~ -0.1. These results 
were mo stly in agreement with a simulation cam paig n con- 
ducted bv lMaccio. Dutton. & van den Bosch I d2008l) and lMaccio et 
(l2007h who covered a mass range 10 9 - 10 13 /z 'M Q , although 
with a slight discrepancy (~ 10%) in the normalization. iDuffv et aL 
(12008b carried out another set of simulations with three different 
box sizes (25, 100 and 400 /T 1 Mpc), each with 512 3 particles 
covering a mass range of 10 11 - lO^/r'M© using the best-fit 
WMAP5 cosmology. They concluded that the median c-M 
relation is lower by about 23% at the low mass end and 16% 
at the high mass end compared to the MS results in the mass 
range of 10 11 - lO^/r'M^. In yet another set of simulatio ns, 



iKlypin. Trujillo-Gomez. & Primackl (1201 Oh and lPrada etaD (1201 lh 
have claimed that the concentration, instead of flattening out at 
high mass, in fact rises. 

Given this context, our primary purpose is to improve the sta- 
tistical power in determining the c-M relation and its scatter at 
high masses, while retaining good mass resolution, and second, 
to study the behavior as a function of redshift and cosmology. 
Finally, we note that the improved statistical power is important 
in comparing with observations of massive clusters as the num- 
bers of well-observed clusters is expected to rise significantly 
in the near future (in the past, simulations may have contained 
only one cluster at the upper mass end, where we have hun- 
dreds). 

3. SIMULATION SUITE 

Throughout this paper, we use the following ACDM cosmol- 
ogy as a reference: u> m = 0.1296 (Ct m = 0.25), uib = 0.0224 {ilt, = 
0.043), n s = 0.97, cr 8 = 0.8, and h = 0.72 where u = VLh 2 and 
f2 m represents the total (dark + baryon) matter density. We as- 
sume spatial flatness. This model is in excellent agreement with 
the latest best- fit cosmological mod el provided by WMAP-7 
measurements dKomatsu et al. lIToi ll) . In order to cover a wide 
range of masses, we analyze three simulations with different 
volumes and number of particles. A summary of the runs is 
given in Table Q] The mass resolution in the large-box run 



is sufficient for measuring the concentrations for halo masses 
> 10 14 /T 1 M Q , with a minimum of 2000 particles per halo. At 
z = 0, we have more than 100,000 such halos, therefore our 
statistical control may be considered to be more than satisfac- 
tory. In the MS and lDuffv et al~l J2008I) simulations, the largest 
boxes used are of size 500/j _1 Mpc and 400/r'Mpc respectively, 
with limited statistics for cluster size halos in the mass range 
10 14 - 10 15 /i _l Mo. We provide a large sample of cluster size ha- 
los, with roughly 64 times more volume th an in the MS run an d 
1 25 times more than in the simulations bv lDuffvetaTI (120081) . 

The largest simulation (both with respect to volume and par- 
ticle number) is carried out using our new Hardware/Hybrid 
A ccelerated Cosmolo gy Code (HACC) fram ework described 
in lHabib et all (120091) and lPope et al. I (120101) . This simulation 
covers a volume of (2 /2~'Gpc) 3 and evolves 2048 3 particles and 
was run on the hybrid supercomputer Cerrillos at Los Alamos 
National Laboratory. (Another 2048 3 particle run with a 512 
/r'Mpc box was used to test the results obtained at z = 2 from 
the GADGET-2 run described below.) The HACC framework 
has been designed with flexibility as a prime requirement; it is 
meant to be easily portable between high-performance comput- 
ing platforms based on different architectures. The first version 
of the code has been optimized to run on the Cell-hybrid archi- 
tecture shared by Roadrunner (the first computer to break the 
Petaflop barrier) and Cerrillos. A first extension of this version 
of the code has been developed for hybrid CPU/GPU systems, 
written in OpenCL. 
al. I HACC's code structure is split into two components: a long- 
range force solver and a short range module. The long-range 
force solver uses a parallel Particle-Mesh (PM) alg orithm with 
spect ral filtering and super-Lanczos differentiation jHamming I 
Il998l) . In this part of the code, the long-range force is calcu- 
lated by depositing tracer particles onto a regular grid and us- 
ing Fourier transform methods to solve the Poisson equation 
(with in effect a modified Green function) and then interpolat- 
ing the force from the grid back onto the particles. The spectral 
component of the code is implemented in C++/MPI and can 
run on any standard parallel machine. The current 2-D domain- 
decomposed implementation of the FFT allows it scale to mil- 
lions of MPI ranks. The implementation of the particle depo- 
sition and force interpolation routines depends on the machine 
architecture. On Roadrunner and Cerrillos, these routines were 
implemented on the Cell processor. 

The short-range module adds the high-resolution force be- 
tween particles and can be implemented in different ways and 
on different platforms. On Cell and GPU -based systems, anN 2 - 
algorithm is used to evaluate the short range forces (in chaining 
mesh patches), leading to a P 3 M implementation. This works 



Table l 

Description of the simulation suite 



Code 


Box 


Softening 


Particles 


m p 




[/t-'Mpc] 


[/T'kpc] 




[/r 1 


M Q ] 


HACC (HACC) 


2000 


7 


2048 3 


6.5 


10 1U 


HACC (HS) 


512 


7 


2048 3 


1.1 


■10 9 


GADGET-2 (G) 


936 


36 


1024 3 


5.3 


10 10 


GADGET-2 (GS) 


128 


10 


512 3 


1.1 


•10 9 



Note. — Runs are referred to in the paper by the names in parantheses. HS 
was run up to z = 2. 
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well on hardware-accelerated machines since it is computation- 
ally intensive and uses a simple data structure. The P 3 M ver- 
sion of HACC has been e xtensively tested agains t the code 
comparison suite results of iHeitmann et al.l d2005l) . On non- 
heterogeneous systems, such as the IBM BG/Q, the A^ 2 algo- 
rithm is replaced by a recursive coordinate bisection (RCB) tree 
method to guarantee good performance. 

In addition to the main code base, a parallel analysis frame- 
work for HACC has been developed. This framework runs on 
conventional supercomputing hardware (or on the 'top' layer 
of a heterogeneous system). Among other utilities, it contains 
a halo and sub-halo finde r. The halo finder was part of a re- 
cent comparison project dKnebe et al. 1 1201 II) and is used for 
the analysis results presented in this paper. Major parts of the 
HACC analysis framework have been implemente d into Par- 
aView and publicly released (IWoodring et al. ||201 lb . 

The smaller simulations are c arried out with G ADGET-2, a 
publicly available TreePM code dSpringel 112005b . Of these, the 
larger simulation - (936 /r'Mpc) 3 volume, 10 24 3 particles - is 
part of the Coyote Universe simul ation suite dHeitmann et al. I 
l2010ll2009l:rLawrence et al. l2010b which spans 38 wCDM cos- 
mologies. This simulation was also used to derive a high-precis ion 
ACDM mass function prediction dBhattacharva et al~l |201 lb . 
(The HACC mass function in the large-volume run presented 
here is in excellent agreement with these results). The smallest 
of the three simulations - (128 /r'Mpc) 3 volume, 512 3 parti- 
cles - serves three purposes: (i) It allows us to probe halos at 
small masses, (ii) it provides large overlap with previous work 
and therefore connects our new results to a mass range that has 
been extensively studied in the past, and (iii) because it is run 
with a completely different code, it provides an excellent check 
on possible code-related systematics (for which we find no ev- 
idence - more details are in the Appendix). 

The initial conditions for all three simulations are generated 
using the Code for Anisotropies in the Microwave Background 
(CAMB 0) and the Zel'dovich approximation at a high start- 
ing redshift, z, ~ 200. Further discussions regarding simulation 
accuracy issues can be found in the Appendix. 

4. ACDM RESULTS 

4.1. c—M relation 

In our simulations, we iden tify halos using a fast parallel 
friends-of-friends (FOF) finder dWoodring et al. 1201 lb with link- 
ing length b = 0.2. (The Appendix contains a discussion of 
how this choice can affect results.) The effects of maj or sub- 
struct ure, relevant for roughly a quarter of the halos (e.g jLukic et al 
l2009b is checked by using morphological cuts mentioned later 
below. Since we are concerned only with the mass profiles, and 
not the dynamical state of the halo, we do not use any velocity 
information (for instance, whether to unbind particles or not). 
Once a halo is found, we define its center via a density maxi- 
mum criteria - the location of the particle with the maximum 
number of neighbors. This definition of the halo center is very 
close to that given by the potential minima. Given a halo cen- 
ter, we grow spheres around it and compute the mass in radial 
bins. Note that even though an FOF finder is used, the actual 
halo mass is defined by a spherical overdensity method, con- 
sistent with w hat is done in observations. (Fo r di scussions on 
halo m ass, see lWhite ll200lllLukic et al~ll2009l and lMore et al~l 
1201 1[ ) Although the mass could be measured independently of 

1 http://camb.info 



the concentration we fit both together to the halo profile, as this 
is potentially less sensitive to fitting bias. (In practice it makes 
little difference.) 

We write the NFW profile as 

C>Pcrit 



P(r) = 



(1) 



(r/r s )(l + r/r s ) 2 

where 8 is a characteristic dimensionless density, and r s is the 
scale radius of the NFW profile. The concentration of a halo 
is defined as ca = ca/ r s , where A is the overdensity with re- 
spect to the critical density of the Universe, pdt = 3H 2 /8irG, 
and ta is the radius at which the enclosed mass, Ma, equals 
the volume of the sphere times the density Ap cr ; t . We compute 
concentrations at two radii corresponding to A = 200 and A = 
A vir , corresponding in turn to c 2 no = Rioo/r s and c vir = R vir /r s . 
The value of A vi > is given by the spherical top-hat collapse 
model; it changes with redshift and cosmology and, for ACDM, 
can be approximated by a fitting formula A v i r = 187r 2 + 82x- 

39x 2 with x = n m (z)- 1, njz) = n m (i+z)7(a„(i+z) 3 + oo 

dBrvan & Norman 1 11998b . For our reference cosmology, A V! > 
varies from ~ 95-170 over the range z = 0-2. We also provide 
a fit for the overdensity of 200 times the mean density, pt, of 
the universe at a particular redshift, z- Written in terms of the 
critical density, this corresponds to A = 200£l m (z) which varies 
from 50-180, over the range z = 0-2 for our reference cosmol- 
ogy- 

The mass enclosed within a radius r for an NFW halo profile 
is given by 

M(<r)=^fl^M A , (2) 
m(CA) 

where m(y) = ln(l +y) — y/(l +y). The mass in a radial bin is 
then 

M 1 =M«r,)-M«r,_ 1 ). (3) 

We fit Eq. [3] to the mass contained in the radial bins of each 
halo, by minimizing the associated value of x 2 as 

(M? m -Mj) 2 



x 2 =E 



(Mf m ) 2 /rn 



(A) 



where the sum is over the radial bins, «, is the number of parti- 
cles in a radial bin, Mf m is the mass in bin i calculated from the 
simulations and M, is the mass calculated assuming the NFW 
profile. The advantage of fitting mass in radial bins rather than 
the density is that the bin center does not have to be specified. 
Note that we explicitly account for the finite number of parti- 
cles in a bin. This leads to a slightly larger error in the profile 
fitting but minimizes any possible bias due to the finite number 
of particles, especially near the halo center. 

We fit for two parameters - the normalization of the pro- 
file and the concentration. Halo profiles are fitted in the ra- 
dial range of approximately (0.1 - l)R V i,- This choice is mo- 
tivated partly by the observations of concentrations that typi- 
cal ly exclude the cen tral region of clusters (e.g., observations 
bv lQguri et aT~ll20Tll to which we compare our results in Sec- 
tion |S]l. More significantly, however, this excludes the cen- 
tral core which is sensitive to the effects of baryonic physics 
and numerical errors arising from limitations in both mass and 
force r esolution, as discuss ed in the Appendix. As already men- 
tioned, [Duffy^etaQdlonJ have shown that, at r < 0AR vir , halo 
profiles are sensitive to the impact of baryons with the profiles 
being affected at r = 0.05/?,,,> by as much as a factor of 2. In the 
Appendix we discuss the robustness of the obtained c-M rela- 
tion as the fitting range is varied; we find that different fitting 
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FIG. 1 . — c-M relations at radii r = R200 and r = R vlr for z = (black), 1 (red), and 2 (blue) for the full (left panels) and relaxed samples (right panels), combining 
results of multiple simulations. The black solid lines at z = are power law fits , a = -0.08 for the full sample and a = -0.084 for the relaxed halos. The solid red 
and blue curves are from the global fit (across all redshifts) discussed in Section l4~2l and shown in Fig. [2] The error bars represent the error in determinin g the mean 
of the concentration in each mass bin. At a given mass, the distribution of concentrations is Gaussian with standard deviation q- c /c~ 1/3 ( Cf. Se ction B31 - the 
shaded region shows the lcr boundary for z = 0. The dotted curves are fitting formulae for the median concentration as given bv lDuffv et a!~l <2008T) . 
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FIG. 2. — c-v relations at radii r = R200 and r = K„iV for z = 0, 1, and 2, for the relaxed and full samples where v = S c /cr(M/^) with 5 ~ 1.676, varying only mildly 
with redshift. The lines are global fits to the data points using a simple assumption for redshift evolution. 



Bhattacharya, Habib, Heitmann, & Vikhlinin 



7 



ranges - chosen with a fair degree of lattitude - agree with each 
other to better than 10% accuracy (Figure lAl U . 

The c-M relation is calculated by weighing the individual 
concentrations by the halo mass, 

EiCiMi 



c(M) = 



(5) 



where the sum is over the number, Nu of the halos in a mass 
bin. The mass of the bin is given by 

M = Y / M i /N i . (6) 

The error on c(M) is the mass-weighted error on the individual 
fits plus the Poisson error due to the finite number of halos in 
an individual bin added in quadrature, 



Ac(M) : 



X/ Ac/My V c 2 (M) 



(7) 



where Ac, is the individual concentration error for each halo. 
The first term dominates towards the lower mass end where the 
individual halos have smaller number of particles and the sec- 
ond term dominates towards the higher mass end, where there 
are fewer halos to average over. 

Figure Q] shows the mean c-M relation obtained from our 
simulation runs weighted by the mass (Eq. [5). We show the 
c-M relation both for relaxed halos and for the full (relaxed + 
non-relaxed) sample. 

To select the relax ed s ample we use criteri a similar to those 
of iNeto et al.1 d2007l) and lDuffv et al~l d2008l) . defining relaxed 
halos as those in which the difference between the location 
of the center-of- mass and the center density maximum is < 
0.07R vir (see also lThomas et al. l200ll) . lNeto et al. Id2007l) have 
used two additional criteria to select their relaxed sample but 
found that the difference in the center of halos method already 
selected most of the relaxed sample. We do not impose their 
additional criteria as it would lead to insignificant changes in 
our sample selection. At z = 0, the relaxed fraction varies from 
0.73-0.6 from M 200 = 10 12 -7.5 x 10 14 /j _1 M q , an d the results 
for thi s fraction are consistent with those found by INeto et al. I 
d2007h . As the redshift increases, one would expect this ra- 
tio to decrease as a function of mass. For the bins centered at 
M 200 = 2.47 x 10 12 /i _1 M Q , andM 20 o = 1.39 x 10 13 /z^Mq, the 
values are 0.77, 0.69, 0.67, and 0.74, 0.63, 0.63, at z = 0, 1, 2 
respectively. AtM 2 oo = 1.39 x 10 14 /j _1 M q , the values are 0.63, 
0.48, at z = 0, 1 (insufficient statistics at z = 2). 

From Fig. Q] it is clear that the c-M relation becomes con- 
siderably flatter at z > 0, with the full sample relation flattening 
more at higher redshift compared to that for the relaxed sample. 
The c-M relation for the relaxed sample has on an average a 
10% higher amplitude compared to that for the full sample. The 
c-M relation at the radius corresponding to A = A n > has about 
a 30% higher amplitude compared to that at A = 200. 

Because the cosmologies considered are essenti ally the same, 
we can directly compare our results with those of iDuffv et al. I 
d2008l although their statistics become somewhat limited near 
the upper end of halo masses. We find that at z = 0, at cluster 
mass scales, their c-M amplitude is about 15% lower com- 
pared to our results. At z = 2, the results from IDuffv et al. I 
(120081) are about 15% higher, but with significant scatter. In 
general, their redshift evolution appears to be slightly com- 
pressed, more so in the case of relaxed halos. Within the sta- 
tistical limitations mentioned, we may consider the compari- 
son to be quite reasonable. At z = 0, our results can be fit- 
ted very accurately by a power law with the exponent, a = 



-0.08 and -0.084 for the full and the relaxed sample. The 
logarithmic sl ope corresponding to the full sample is precisely 
that found by lHavashi & Whiteld2008T) using the halo-density 
cross-correlation applied to data from the MS. The normaliza- 
tion, however, is not expected to be the same because of the high 
value of us = 0.9 chosen for the MS. Note that three different 
analyses of the MS (and an associated smaller-volume, higher 
mass resolution run) have produced slightly discre pant results, 
differing from each other at the 10-20% lev el dNeto et al. I 
l2007l:lGao et al. 120081: lHavashi & White 120081) . This is proba- 
bly a useful empirical measure of the systematic issues inherent 
to halo selection and fitting. In general, the MS results are con- 
sistently higher at all redshifts by about 15% at z = to about 
30% at z = 2 largely because of the higher value of erg. 

4.2. The c — v relation 

We find that the c-M relation becomes almost flat at higher 
redshi ft, with C200 ~ 3 in a way similar to the fi ndings of Gap et al.l 



(12008b (although for a higher erg in their case). IGao et al. Id2008l) 
use the Einasto profile for stacked halos - with its extra shape 
parameter - to account for the flatter c-M relation at higher z. 
Alternatively, we ask if it is possible to explain the flattening 
of the c-M relation with redshift using only the NFW profile 
without adding an extra parameter. To do this, we follow the 
strategy adopted in parameterizing halo mass functions and in- 
vestigate the concentration measurements as a function of the 
rms density fluctuation er(M, z), rather than M (for each of the 
three halo mass definitions). As the central variable, we use the 
peak height parameter, v = 6 c {z) / cr{M, z), where S c (z) is the lin- 
ear collapse threshold. (6 C = 1 .673 for the reference cosmology 
and varies only mildly with cosmology and redshift.) <r(M,z) 
specifies the variance of the matter fluctuations over the scale 
oc M 1 / 3 at a redshift z. As shown in Figure |2j the shape of the 
c—v relation is approximately constant over the redshift range 
z = 0-2, in contrast to the shape of the c-M relation. 

Overall, the evolution of cnooiv) proceeds as ~ D(z)° 5 where 
D(z) is the linear growth factor at redshift z, or by about 30% 
from z = 0-2. The overall z-evolution in our work i s roughly 
consis tent with the z-evolution seen in the MS result of lGao et al. I 
d2008l) . The evolution of c V i r (v) possesses a somewhat larger 
dynamic range and the evolution goes as ~ D{z). The slope of 
the c—v relation is slightly larger for A = 200 compared to that 
for A = A V! >. The amplitude of the c — v relation is only a little 
larger for the relaxed sample compared to the full sample (by 
about - 10%). 

Fitting formulae for c{v) as derived from the simulations for 
the reference ACDM cosmology are given in Table [2] We also 
provide an approximate fitting formula relating v and M valid 
for all overdensities, redshift, and cosmology, which can then 
be used to convert the relations for c - v to those for c-M. 
Table[2]also provides the c—v relation for A = 200S!,„(z) corre- 
sponding to an overdensity of 200p/,. At z = 0, for the reference 
cosmology used here, 200p/, is about half of the virial overden- 
sity. Consequently, the c — v amplitude is about 30% higher 
compared to the amplitude at the virial density. At z = 1 and 
2, the mean density and the virial density become comparable. 
Thus the c—v relation for A = 20051„,(z) has more z-evolution 
when compared to that for the virial density. 

4.3. The distribution of concentrations 

The mean c-M relation needs to be augmented with a good 
quantitative understanding of the concentration scatter around 
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the mean, especially at cluster-scale masses, where simulations 
have historically suffered from lack of volume coverage. Pre- 
dictions for the distribution of the concentration are particularly 
valuable since they can be used to check for selection biases in 
observations. As an example, if there is a concern that lensing- 
based searches are likely to be biased towards high concentra- 
tion halos, then, at a given mass bin, one can test for this bias 
by comparing to the predicted theoretical distribution. For this 
method to work, there should be enough objects at a given mass, 
a target that will be attained in the near future. 

We have computed the concentration distribution for a large 
set of cosmol ogies, a subset of which we discuss here. Pre- 
vious studies (|Jing1l2000b [Shaw et al. l l2006tlNeto etal. Il2007t 
iDuffv et al. 112008b . have fitted the concentration distribution to 
a log-normal distribution. However, this distrib ution is also 
very w ell described by a Gaussian as noted by iLukic et al. I 
d2009b and lReed. Koushiappas. & Gaol d2011l) . We have found 
that a Gaussian distribution provides a very good fit to our data, 
with relatively small non-Gaussian tails. As a representative 
example, we show the distribution of C200 (full halo sample at 
z = 0) for the mass bin centered at 1 .5 x IQ 14 Ii~ 1 Mq in Figure[3] 

Assuming a Gaussian distribution, the standard deviation in 
the c-M relation can be calculated as, 



a c (M) : 



-c(M) 2 



(8) 



and the associated error in determining the scatter is the Poisson 
error in each bin, 



(9) 



where N/ is the number of halos in the mass bin with mass M. 

As illustrated in Figure [3] for the case of halos in the mass 
bin at 1.5 x 10 14 /r'M©, at z = 0, the standard deviation of 
the Gaussian distribution is roughly <r c = 0.33c, over our mass 
and redshift range, with mild dependence on the mass at the 
very high mass end, such that for M200 > 8 x 10 14 /t'Mq, 
<t c ~ 0.28c (Figure|4]i. These resul t s are i n very good agreement 
with iReed. Koushiappas. & Gaol (|201 lb who find a c ~ 0.28c 
for an analysis of halos extracted from the MS. If instead, for 
comparison purposes, we fit our concentration distribution us- 
ing the log-normal function, we get cr(log 10 (c2oo)) = 0-16 for 
the full sample and 0.12 for the relaxed sam ple. This may be 
compared to the result of IDuffv et al. I ( 20081). w ho obtain 0.15 
and 0.11, respectively, and to lNeto et al. I d2007b who find 0.14 



Table 2 

c{y) — U FITTING FORMULAE. 





A = 200 


A = A Vi> 


A 


= 200,% 


full 


£>(z)°- 54 5.9^ a35 


D(z) a9 7.7^ - 29 


D(z) 1 


.15 9 ^.29 


relaxed 


D(z)°- 53 6.6i/- a41 


D(z) 101 8.9i/- a34 


D(z) 1 


2 io.i^-° 34 



Std. Dev. 



v-M 



a c = 0.33ca 



v(M,z) « ggj 



and 0.1. Ou r scatter is therefor e ab out 5 - 10% higher than 
the results of IDuffv et al. I d2008b and iNeto et al. I d2007b . As 



shown by lNetoetal.ld2007h . the variance is at a minimum for 
the radial range (0.05- 1 )/?,,,>. As noted earlier, we fit the halo 
profile over the range (0.1 - l)i?v; r which may account for the 
~ 7% larger value of the standard deviation. Our choice trades 
off a slight increase in scatter for robustness against system- 
atic effects, as previously discussed. Figure [4] shows that the 
relation <j c = 0.33c is more or less independent of mass, red- 
shift, or the dynamical state of the halo. The relation also re- 
mains constant when the cosmology is varied (Figure [6]). This 
means that the standard deviation of the concentration distribu- 
tion depends on cosmology, redshift, or the dynamical state, in 
the same way as the mean co ncentration, confirming the initial 
finding o flDolag et all d2004b from a small sample of simulated 
halos, but spanning multiple dark energy cosmologies. 

5. C-M RELATION FOR wCDM COSMOLOGIES 

In this section we study how the halo profiles, and hence the 
concentration, vary with cosmology. We use 18 different runs, 
each with a volume of ( 1 .3 Gpc) 3 and 1024 3 particles. The runs 
are carried out using GADGET-2 and each run has a differ- 
ent set of wCDM parameters. The se simulations are a su bset 
of the Coyote Univer se suite (see iHeitmann et al. 1 120091 and 
iHeitmann et al. Il2010l for details). The simulation suite con- 
sists of 38 runs covering the 2a range of w CDM parameter 
space as constrained by WMAP 5 year results dKomatsu et al. I 
12009b . We choose 18 runs (plus the reference cosmology run) 
out of the 38 to show the cosmology dependence and retain 
the model numbering from the original Coyote runs. The runs 
have a coarser force resolution than the GS and HACC simula- 
tions. The effect of this is considered in the Appendix, where 
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FIG. 3. — c-M distribution at a mass bin centered at 1.5 X 10 14 h~ Mq using 
results from the HACC simulation at z - 0. Lines show a Gaussian distribution 
with standard deviation a c /c ~ 0.33. 
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Table 3 

Parameters for the 18 cosmological models used to study the c-M relation 



# 




u m 




n s 


— w 


ff8 


h 


M 

10 14 M Q 


variation 

% 


1 





1539 


0.0231 


0.9468 


0.816 


0.8161 


0.5977 


13.3 





3 





1324 


0.0235 


0.9984 


0.874 


0.8484 


0.6763 


9.96 


+15 


4 





1381 


0.0227 


0.9339 


1.087 


0.7000 


0.7204 


4.42 


-18 


5 





1358 


0.0216 


0.9726 


1.242 


0.8226 


0.7669 


7.20 


-20 


7 





1268 


0.0223 


0.9210 


0.700 


0.7474 


0.6189 


7.30 


+15 


8 





1448 


0.0223 


0.9855 


1.203 


0.8090 


0.7218 


8.04 


-15 


9 





1392 


0.0234 


0.9790 


0.739 


0.6692 


0.6127 


4.98 


+10 


12 





1223 


0.0225 


1.0048 


0.971 


0.6271 


0.7396 


2.26 


+14 


13 





1482 


0.0221 


0.9597 


0.855 


0.6508 


0.6107 


4.78 


+5 


1 J. 
1 -+ 


n 
u 


1/L71 
1*+ / 1 




1 rnnfi 

1 .WjuO 


1 .U1U 


O 707 s 
\J. 1 \J 1 J 


W.DDoo 




n 
u 


15 





1415 


0.0230 


1.0177 


1.281 


0.7692 


0.7737 


5.47 


-15 


16 





1245 


0.0218 


0.9403 


1.145 


0.7437 


0.7929 


4.22 


-10 


17 





1426 


0.0215 


0.9274 


0.893 


0.6865 


0.6305 


5.50 





18 





1313 


0.0216 


0.8887 


1.029 


0.6440 


0.7136 


3.05 


-10 


19 





1279 


0.0232 


0.8629 


1.184 


0.6159 


0.8120 


1.88 


-15 


20 





1290 


0.0220 


1.0242 


0.797 


0.7972 


0.6442 


8.24 


+20 


30 





1234 


0.0230 


0.8758 


0.777 


0.6739 


0.6626 


4.09 


+10 


37 





1495 


0.0228 


1.0233 


1.294 


0.9000 


0.7313 


11.7 


-15 



Note. — The second column from the right shows the mass corresponding to v = 3 for each cosmology at l = 0. The right-most column shows the approximate 
variation of the mean c—M relation with respect to the reference ran. 




FIG. 4. — The t'2oo ~^200 distribution for the relaxed and full halo samples 
characterized by the ratio of the standard deviation to the mean value of C2oo- 
All three redshifts are plotted. Note that cr c /c shows no redshift evolution. The 
case of A = A v ,y shows identical behavior. 



FIG. 6. — The C200 _ ^200 distribution at z = (full sample) when »'CDM 
parameters are varied, following the characterization of Fig. [4] The scatter is 
larger at high masses due to lower numbers of halos in high-mass bins. 
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FIG. 5. — Mean c-M relation at z = (full sample) when vvCDM parameters are varied. The solid curve in all three panels is the fit to the reference ACDM 
cosmology as specified in Table ff] 



it is shown that the G run has a systematically lower concentra- 
tion compared to the HACC run over the same mass scale by 
about 5-10% (Figure lAl II left panel). To compensate for this 
minor underestimate, we rescale concentrations obtained from 
the wCDM runs by a factor of 1 .05, checking for correctness by 
comparing against the fit obtained for the reference cosmology. 

Figure [5] shows the variation of the c- v relation with respect 
to the best-fit WMAP5 cosmology. The mean c — v relation 
varies by about ± 20% over the currently allowed vvCDM cos- 
mological parameter range. Note that v already accounts for 
some of the cosmology dependence of the c-M relation, so a 
part of the variation is actually hidden. Since we have already 
found that expressing c as a function of v explains the redshift 
evolution of NFW halo profiles, we illustrate the cosmology de- 
pendence using the c-v relation in place of the c-M relation. 

Table [3] shows the approximate difference between the (cor- 
rected) mean c-M relation seen in each of the wCDM runs 
compared to the mean c-M relation obtained for the reference 
ACDM cosmology. Note that although most of the variation 
in the c-M relation is in the overall amplitude, the slope also 
changes for some of the models (e.g., M003, M012). Inter- 
estingly, we find that some of the models show no variation 
compared to the reference, although these models differ across 
the range of cosmological parameters. For example, M014 and 
MO 17 both have lower 0% compared to the reference model, but 
show essentially no variation - parameters other than erg are 
clearly also active. The standard deviation of the concentra- 
tion distribution, on the other hand, changes in the same way as 
the m ean, leaving the ratio a c jc almost universal (IDolag et al. I 
12004 . Figure [6] shows that the a c /c varies by < 5% over the 
range of vvCDM cosmologies. 

Sem ianalytical 'toy models' b ased on Press-Schechter argu- 
ments dPress & Schechtenll974l) have attempted to model the 
cosmology and redshift dependence of halo concentrations via 
the underlying dependence on the matter power spectrum and 
the evolution history of the universe (iNavarro. Frenk. & White 



[l997t 



Eke et al. l200ltlBullock et alj2001l) . The model o flNavarro. 



dl997l) has two free parameters. It defines the halo formation 
redshift by requiring that half of the final halo mass be in pro- 
genitors with masses of some fraction of the final mass. The 
characteristic density scale of the NFW profile is then set by 



assuming it to be some multiple of the cosmic density at the 
redshift of halo formation. The mass fraction a nd the density 
multip liers are given by fitting to simulations. iBullock et al.l 
j2001l) modified this prescription by redefining the formation 
redshift as the redshift where the nonlinear mass scale M* is 
some fraction of the final halo mass. They predicted the concen- 
tration as a multiple of the ratio of scale factors at the formation 
and collapse redshifts. Again, the fraction and mu ltiplier are 
floatin g parameters, determined by fitting. Finally, lEke et akl 
d2001l) used a single parameter (calibrated using simulations) 
to connect the collapse redshift with the effective amplitude of 
the power spectrum at cluster scales. 

Even though the models reproduce the expected behavior 
of the c-M relation discussed in the Introduction, quantita- 
tively their success has been decidedly mixed - results have 
been satisfactory over limited dynamic ranges when fitted to 
simulations, but tended to break down as the range is extended 
or cosmological parameters are varied. Particularly significant 
for us is the breakdown at large halo masse s , characteristic of 
massive clusters at z = dNeto et al. Il2007t iGao et al71 l2Q08h 
iDuffv et al. l2008tlziia"o et al. 120091) . with a corresponding break- 
down at lower masses, but at higher redshifts. Additionally, the 
redshift evolution of the c-M relation in the models is signifi- 
cantly stronger than is actually seen in simulations. 

Althou gh we have not a ttempted to optimize model param- 
eters, the Eke et al.l d2001l) prescription agrees with our simu- 
lation results at 10-20% accuracy for C200 for the ACDM cos- 
m ology (see Table [2] f or the fit). We find that the prescriptions 
of IDolag et al. I d2004l) (using a growth factor ratio multilpier) 
do not explain the cosmology variation in concentrations that 
we observe. For instance, in our case, the growth factor only 
varies by < 5% over the range of simulations used, whereas the 
variations of concentrations seen is ~ ± 20%. 

Regarding the modeling of the scatter in the concentration, it 
is natural to examine this in the contex t of different assembly 
histories for halos with the same mass dWechsler et al. 112002b 
Fte«Bo&t\Mhllad03l) ( See also. ICohn & White 1120051 ) However, 
in their MS analysis. iNeto et al7ld2007l) find that the concentra- 
tion scatter cannot be accounted for by differences in the time 
of form ation alone. Additiona l consequences of environmental 
effects dWechsler et al. I [2006b appear to be important primar- 
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ily for low-mass halos. Therefore one is driven to the general 
conclusion that there is still no replacement for large-scale sim- 
ulations if reliable predictions for halo concentrations and the 
distribution of concentrations are required. 

6. COMPARISON WITH OBSERVATIONS 

In this section, we compare our simulation results with some 
of the recent observations of the concentration-mass relation for 
clusters. The observational results span a variety of techniques, 
including strong and weak lensing (e.g., Comerford & Nataraian I 



( 2007| have measured the concentration of 34 dynamically re- 
laxed clusters (0.06 < z < 0.69) from Chandra observations 
(left panel). The theoretical predictions are in good agreement 
for masses M,,,> > 4 x 10 14 /r' Mfn, with minor ten sion at lower 
masses. The data presented by lBuote et al. I (l2007h are a compi- 
lation of analyses of relaxed systems from Chandra and XMM- 
Newton; we sho w only the higher mass range, rep resent ed by 



2007]: [B: 



Broadhurst et alj2008l:lMandelbaum et al. l2008HOkabe et al. 
Qguri et aI71 201 ll) X-ray observations of relaxed clusters 



(e.g.. Buote et al. l2Q07hl Vikhlinin et al. l2006l:ISchmidt & Allen I 



120071: 1 Vikhlinin et al. Il2009h and relaxed and unrelaxed clusters 
dEttori et al. 1201 1[). and cluste r kinematics (e.g jRines & Diaferiol 
l2006t IWojtak & Lokas 1 1201 Oh . Our aim is to provide a set of 
figures that enables the reader to judge by eye the current status 
of how well the theoretical predictions match against observa- 
tions. Because there are significant observational systematics 
that are unclear and the observational statistics are still limited, 
we do not believe that a more complete statistical analysis is 
necessary, or even particularly useful. The strategy we follow 
is to take the ratio of each measured concentration to the theo- 
retically predicted concentration at the object's observed mass 
and redshift. We then bin in mass to show a relatively limited 
number of comparison points in each figure. Thus each point in 
the observation plots represents an average over ^5 data points. 
(The corresponding Poisson error bars use the improved for- 
mula a± = ^Nh + l/4± 1/2 as given bv lHeinrichll2003L which 
asymptotes to y/NJ, at large Ni,.) 

We begin our comparison using results fr om X-ray observa- 
tions of relaxed clusters as shown in Fig. [7] I Schmidt & Allen I 
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FIG. 8. — Ratio of observed concen tration to theoretical predictions for the 
XMM-Newton cluster observations of Ettori et al. 1 1201 II) . The shaded area 
represents the la boundary for the theoretical predictions. Each data point 
actually represents observations of multiple clusters. 



results taken fro m Pointecouteau, Arnaud, & Pratt (2005) . | Vikhlinin et al 
d2006l) . IZappacosta et al. Id2006l) and lGastaldello et al. I (120071) . 
spanning a redshift range of 0.016 < z < 0.23. All of these 
] results are in very good agreement with the predictions, lying 
within the 1<t boundary. The right panel shows results from 
19 clust ers that were part of t he Chandra Cluster Cosmology 
Project dVikhlinin et al. 112009b (0.016 < z < 0.25), the dataset 
represented in Table [5] Once again, the agreement is excellent. 
Overall, we conclude that comparisons with X-ray measure- 
ments of relaxed clusters are in good accord with (concordance) 
ACDM predictions. 

Next we turn to the results o f lEttori et al. Id2011l) who mea- 
sured the concentrations of 44 X-ray luminous clusters (0.09 < 
z < 0.31) using XMM-Newton (Fig. [8}. Their sample contains 
both re laxed and unrelaxed clusters. As with the I Schmidt & Allen] 
d2007l) comparison, we find that the simulation results are in 
good agreement with these observations forA/200 > 4 x 10 14 /i _1 M Q . 
As the authors themselves note, a slope cannot be fitted to their 
data because of the narrow mass range of the observations rel- 
ative to their errors. Thus, we regard the current level of agree- 
ment as being quite satisfactory. 

We now consider lensing measurements of cluster profiles 
using weak and strong lensing and combinations thereof. Fig- 
ure|9]shows the comparison of the theoretical predictions against 
the results of LocUss, a weak lensin g study of 30 cluste rs with 
Subaru/Suprime-Cam imaging data dOkabe et al. WE) 

and a 

combined strong and weak l ensing analysis of 2 8 clusters from 
the Sloan Giant Arcs Survey dOguri et al. 1201 ll) . The left panel 
of Fig. [9] shows the weak lensing results displaye d in the same 
manne r as for the X-ray datasets. The results from lOkabe et al. I 
d2010h are in excellent agreement with our predictions, com- 
pletely consistent with the corre sponding me a surem ents from 
relaxed clusters. The results of lOguri et al. I d201 lh are con- 
sistent with our predictions for M vir > 4 x 10 14 /t'Mq, but at 
lower masses, there appears to be a significant discrepancy, with 
a much steeper c-M dependence. Although baryon cooling 
may play a role at smaller masses, there is no convincing rea- 
son for such a large effect - for whic h there is no signa l in the 
X-ray data (nor in the simulations of iDuffv et al. Il2010h . Note 
th at the target selection in the two surveys is quite different, that 
of lOkabe etaTld2010l) being essentially volume-li mited, while 
any str ong-lensing selected sample such as that of lOguri et al. I 
d201 ll) m ust have a significant amount of selection and projec- 
tion bias dRozo et al. "ll2008UMeneghetti et al. l2010h . Note also 
that an analysis based on mock weak lensing observations in 
the MS dBahe et al. 11201 ll) has shown that there is bias in weak 
lensing measurements of concentration as well, tending to de- 
press the measured concentration by a small amount from the 
predicted value. 

The right panel of Fig. [9] shows the combined stro ng plus 

weak lensing analysis including a model for lensing bias dOguri et al. I 
l20Tll) . Processing the results through the lensing bias model 
(by enhancing the theoretical prediction) brings down the dis- 
crepancy significantly but there is still evident tension for masses 
Mvir < 4 x 10 14 /2 _1 M(7,. Nevertheless, we note that there is a 
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Table 4 



OBSERVATION DATA USED IN THIS PAPER 



observation 


method 


rel./all 


# clusters 


redshift range 


Oguri et. al. 


Strong+Weak lensing 


all 


28 


0.28 < z < 0.64 


Okabe et. al. 


Weak lensing 


all 


30 


0.15 <z< 0.3 


Wojtak & Lokas 


Kinematics 


rel. 


41 


z<0.1 


Vikhlinin et al 


X-ray 


rel. 


19 


z<0.2 


Schmidt & Allen 


X-ray 


rel. 


34 


0.06 < z < 0.7 


Buote et. al. 


X-ray 


rel. 


26 


z < 0.23 


Ettori el. al. 


X-ray 


all 


44 


0.1<z<0.3 



Note. — We use only objects with mass > 5 X 10 13 h 'Mq from Buote et al. 



Table 5 

Updated masses and concentrations from the Chandra Cluster Cosmology Proiect 



Cluster 


M 500 (M ) 


SM (M ) 


C500 


+6 C 




6c 




z 


al33 


3 


166xl0 14 


3 


776 xlO 13 


3 


15 





29 





28 





0569 


a262 


8 


310x10° 


7 


272xl0 12 


3 


48 





30 





30 





0162 


a383 


3 


049 xlO 14 


3 


100xl0 13 


4 


31 





42 





40 





1883 


a478 


7 


668xl0 14 


1 


OlOxlO 14 


3 


57 





27 





26 





0881 


a907 


4 


623xl0 14 


3 


790xl0 13 


3 


46 





42 





42 





1603 


al413 


7 


569xl0 14 


7 


550xl0 13 


2 


93 





18 





17 





1429 


al795 


6 


009 xlO 14 


5 


134xl0 13 


3 


21 





18 





18 





0622 


al991 


1 


235 xlO 14 


1 


654xl0 13 


4 


31 





34 





34 





0592 


a2029 


8 


147xl0 14 


7 


674xl0 13 


4 


04 





21 





21 





0779 


a2390 


1 


077xl0 14 


1 


092xl0 14 


1 


66 





13 





13 





2302 


clll59 


1 


056xl0 14 


2 


578xl0 13 


1 


77 





38 





24 





0810 


MKW4 


7 


734x10° 


1 


032xl0 13 


2 


54 





16 





14 





0199 


a2717 


1 


478xl0 14 


2 


134xl0 13 


2 


69 





19 





19 





0498 


a3112 


3 


448 xlO 14 


3 


097 xlO 13 


4 


47 





28 





27 





0761 


al835 


1 


245 xlO 15 


1 


342xl0 14 


2 


81 





17 





17 





2520 


a 1650 


4 


683xl0 14 


1 


736xl0 13 


3 


74 





19 





19 





0846 


a2107 


2 


361xl0 14 


3 


928 xlO 13 


3 


38 





28 





25 





0418 


a4059 


3 


496xl0 14 


2 


691 xlO 13 


2 


95 





09 





09 





0491 


rxjl504 


1 


068xl0 15 


1 


768xl0 14 


3 


16 





38 





38 





2169 



Note. — <5M is the estimated error in the mass, +8 C and — 8 C , are the upper and lower error bounds for the concentrations. The masses are for h — 0.72. 
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clear trend of lensing concentrations reducing over time and be- 
coming more consistent with the theoretical predictions. Other 
data our results appear to be in agreement with can be found in 
IComerford & Nataraian I d2007l) (strong lensing) and lCoe et al. I 
d2012l) (strong and weak lensing). A cautionary note regarding 
weak lensing conce ntration measurements of clust ers is pro- 
vided in Figure 6 of IComerford & Nataraian ( 20071) regard ing 
Abell 1689 and in the results given in Israel et aTl d201 lb as 
part of the 400d weak lensing survey. 

Instead of using individual objects, a stacked statistical anal- 
ysis can be applied to clusters, as carried out usi ng the Sloan 
Digital Sky Survey by iMandelbaum et al. I d2008l) . to a cluster 
mass range of ~ 6 x 10 14 /i _1 M Q . This analysis sees no evi- 
dence for a major boost in concentration at lower masses and 
the final result - C2006 ~ 4.6 ± 0.7 at (z) = 0.22 at a mass of 
M 2 oob ~ 10 14 /2 -1 M is 20-40% less than our prediction of 
cioob ~ 6.5 at the corresponding mass. The mild c-M depen- 
dence they observe is however in good agreement with our pre- 
dictions - ~ 0.09 compared to the observed slope of 0.13 ± 
0.07. 

Finally, we consider the estimate s of the concentrat i on us- 
ing galaxy kinematics in clusters. iRines & Diaferio I J2006h 
matched X-ray cluster catalogs with SDSS and used infall pat- 
terns to compute cluster mass profiles. The C200 concentration 
has significant scatter - values ranging from 5- 17 - but their 
best-fit average profile, with fits restricted to r < R20Q, yields 
C200 = 5.2 which, at an average mass of M200 — 10 14 /r'M^, is 



in agre ement with our predictions. More recent! v lWoitak & Lokasl 
d2010l) have published an analysis of 41 relaxed galaxy clusters 
(0.013 < z < 0.095); we compare our predictions with their re- 
sults in Fig. [TO] Although their data has considerable scatter it is 
in quite reasonable agreement with the predictions from simu- 
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FIG. 10. — Ratio of measured to predicted concentra tions; the data is 
taken from the results of a projected phase space analysis by Woitak & Lokas 
(2010). 



lations. Thus, desp ite possible systematic diffic ulties with such 
methods (see, e.g JWhite. Cohn. & SmitlfeOlOl) . the current re- 
sults are in reasonable accord with theoretical expectations. 

7. SUMMARY AND DISCUSSION 

We presented results for the concentrations of dark matter 
halos using a set of large-volume simulations. With a total vol- 
ume roughly 1-2 orders of magnitude larger compared to pre- 
vious simulations, we focused on studying the c-M relation 
for massive clusters. As shown in the past, at the high mass 
end, the c-M relation becomes flatter at z = and the flat- 
tening becomes more significant at higher z- The mean con- 
centration of the sample when expressed in terms of the peak 
height parameter, v{M,z) = S c /a(M,z), shows a roughly uni- 
form slope at all redshifts. Indeed, the slope of the c-v rela- 
tion does not change with redshift. The amplitude of the c-v 
relation evolves by about 30% at the high mass en d from z = 
0-2. T he z-evolution is consistent with the results o flGao et al. I 
(2008), although the overall amplitude of the concentration dif- 
fers because of the different choices of erg. We do not ob- 
serve a rise in concentration at higher ma sses as reported by 
iKlvpin, Truiillo-Gomez, & Primack l d2010h and lPrada etakl d201ll) 
(the Appendix includes further discussion). Because of our 
large halo sample, we can study the distribution of the con- 
centration in individual mass bins; we find that the distribu- 
tion of concentrations is well-described by a Gaussi an PDF 
dLukic et al. 1 120091: iReed. Koushiappas. & Gao 1 1201 lb . Thus 
the halo profile shape can be described by two parameters - 
the mean concentration and its standard deviation. By compar- 
ing results across a number of wCDM cosmologies, we find that 
the standard deviation is roughly universal, a c = 0.33c, and does 
not change with redshift, halo dynamical state, or cosmological 
parameters. 

We investigated how the concentration changes as the cos- 
mological parameters are varied using a set of 18 runs spanning 
the vvCDM cosmology parameter space. The parameter range 
covers the 2er variation around the best fit WMAP5 cosmology. 
We find that over this parameter range, the c-M relation varies 
by ~ ± 20%, although the standard deviation a c follows the 
relation a c = 0.33c. As suggested by our work on the wCDM 
models, and also previous studies of redshift evolution, the halo 
formation epoch, and hence the concentration, depends on the 
matter fluctuations, slope of the power spectrum and the growth 
factor. Thus calibrating the c-M relation as a function of cos- 
mology is important for a wide variety of problems, ranging 
from galaxy formation, the weak lensing shear power spectrum, 
to the case of assembly bias in clusters. We will address the cos- 
mology dependence of the c-M relation in detail using more 
simulations and analytical models elsewhere. 

The simulation predictions are in good agreement with ob- 
servations from strong lensing, weak lensing, galaxy kinemat- 
ics, and X-ray data for massive clusters with masses M,,„- > 
4 x 10 14 /t'Mq. At lower masses, different observations suf- 
fer from different sources of systematic error. For example, 
the lensing data need to account for bias due to the triaxiality 
of halos while the X-ray data typically ignore the non-thermal 
pressure component in galaxy clusters w hich can lead to a sys- 
tematic underestimate of the cluster mass (lLau et al .120091) . The 
simulations, on the other hand also need to account for baryonic 
effects which play a bigger role as the halo mass decreases. As 
a result, due to cooling, gravity-only simulations may predict 
20 - 30% lower concentration for clusters with masses M,,„- < 
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4 x 10 14 h~ l M Q . The fact that most recent observations are 
in agreement with the simulation results (and amongst them- 
selves) to better than 20% for massive clusters (M VI > > 4 x 
10 14 /2 _i Mq) indicates that baryonic effects influencing the clus- 
ter mass profile are indeed small and that the individual obser- 
vational systematics are under some level of control. 

APPENDIX 

SYSTEMATICS AND ROBUSTNESS CHECKS 

In this Appendix, we investigate various sources of possi- 
ble systematic errors in determining halo concentrations. We 
note that the simulations were carried out keeping in mind er- 
ror control requirements on force-resolution, time-stepping er- 
rors, an d mass resolution that have been spelled out in the lit- 
erature (iTormen Bouchet. & White Ill997t iPower et al. Il2003t 
iHeitmann et al. ll2008tlLukic et al. 112007^ 

We locate halo centers by using a very fast FOF method. In 
principle, there is a very mild selection effect induced by the 
choice of linking length - if the FOF finder links two halos, 
then only the higher density center of one halo will be used, 
an d the other halo w ill be statistically lost (see, e.g. Fig. 12 
in lLukic et al. 112009b . A smaller linking length would in effect 
free up the other halo as well, although it would slow down 
the center-finding algorithm. In practice, however, this is not 
an issue given the fact that individual concentrations need to 
be extracted with small errors. To do so we use a minimum 
number of ~ 2000 particles per halo; tests have shown that with 
more than 500 particles per halo, there is excellent agreement 
between our method of h alo finding and con ventional spherical 
overdensity halo finders (iKnebe et al. 11201 ll) . 

Another problem with halo center finding is miscentering, 
which in general will tend to reduce the concentration. To pro- 
duce a quantitative estimate for this effect, we consider an ex- 
treme example by offsetting the center of every halo by 30/? _1 kpc 
(a value of the same order as the force resolution) and recom- 
puting the profile of every halo in the G run. As shown in 
Figure IA1 II left panel, randomly offseting the center by this 
amount reduces the c-M relation by about 10% at M200 = 5 x 
10 13 /! _1 M Q with the difference vanishing towards the high mass 
end. (This sort of misestimation is more relevant to analyses 
with stacked halos.) 

For profile tests across the different runs, we begin by com- 
paring the GS and HACC runs (left panel of Figure lATTb . Note 
that although these runs have been carried out with completely 
different codes, the data for C200 g° es over smoothly from one 
mass range to another. (GADGET-2 and HACC were run with 
roughly similar force resolution.) The G runs were originally 
carried out for a different purpose, hence their force resolution 
is somewhat lower. As expected, this has the consequence of 
mildly reducing the conc entration dTormen. Bouchet. & White] 
ll997tlPoweret al. 112003b by about 5%. 

We have tested our profile-fitting method by generating Monte 
Carlo NFW profile samples using different numbers of parti- 
cles; with the particle numbers used to sample halos kept larger 
than 2000, the method was accurate to a few percent (worst 
case) and superior to simpler methods based on radius ratios 
and variants thereof. To investigate how other parameters could 
affect concentration values, we went back to using the halos 
from the simulations. The right panel of Figure I A 1 1 1 shows the 
effect of varying the range of the halo profile used to fit to the 
NFW form. Changing the starting radius from r = 0.1/? n > to 
r = 0, with the outer limit fixed at r = R v i r reduces the overall 



c-M relation by about 5% (resolution/particle undersampling 
limitations). Fixing the starting radius at 0.1/? v ; r but changing 
the stopping radius to 2i? V! >, only changes the relation by a neg- 
ligible amount from the fiducial range of (0.1 - l)R V j r . Reduc- 
ing the stopping radius to 0.5J? V! > steepens the c-M relation by 
about 10%. 

Because our primary interest is in halos that have mass sig- 
nificantly in excess of M*, it is important to ask what possible 
systematic effects could arise from fitting such objects with- 
out paying attention to their infall structure. The average ra- 
dial velocity of a halo fluctuates around zero out to an infall ra- 
dius, ri n f, beyond which it goes negative, this transition roughly 
defining the boundary of the infall region. Purely as an informal 
nomenclature, we refer to the region internal to the infall radius 
as the virialized region. We find all halos that have R v , r > r,„/ 
and exclude them from the analysis, thus focusing attention on 
halos that have much less infall contamination. The result is 
shown by the solid line in the right panel of Figure lATTl Not 
unexpectedly, cluster size halos with masses ~ 1O 14 /2 _1 M0 and 
greater are much more sensitive to this cut, and display an en- 
hancement of the c-M relation by about 10% when only the 
'virialized' sub-sample is used. 

As mentioned in Section [2] the halo concentration can be 
measured in different ways, even if one sticks to the NFW defi- 
nition^) of concentration. Therefore, it is important to investi- 
gate what sources of uncertainty can arise from using different 
definitions that may be mathematically equivalent, but not oper- 
ationally the same. Here we study two alternative independent 
techniques for measuring concentrations - the radius ratio and 
the maximum circular velocity. The radius ratio method is very 
simple: We measure the radius at A = 300 and 200 for each 
halo in our sample. Then, assuming the halos are described by 
an NFW profile, the radius ratio is given by 

m(C20oR 3<X)//?200) 



300^oo 
2004 )0 



(Al) 



m{c2oo) 

where m(x) = ln(l +x) — x/(l +x). Given a measurement of 
^30o/^200> one solves the nonlinear equation Eq. lAll for C2oo- 
The mean c-M relation is then obtained in the same way as 
for profile fitting. The left panel of Figure lAT2l shows the con- 
centration measured using the profile fit and by the radius ratios. 
For this test we focused on the relaxed sample, although the full 
sample gave identical results. The diagonal line in the figure de- 
notes the ideal exact agreement of the concentration measures 
from the profile fit and from the radius ratio. As the results in 
the figure show, the two methods agree quite well within the 
specified errors. The error in the concentration measurement 
using the radius ratio arises from the error in determining the 
radii - the Poisson error from the total number of particles in- 
side Z?3oo and /?200- As expected, the radii are determined quite 
accurately as there are large number of particles inside these 
radii, but because of the logrithmic nature of Eq. I A 1 1 the AC200 
become non-negligible. The error on the mean concentration 
also includes the Poisson error due to the finite number of halos 
in the individual mass bins. 

The second method we investigate rel ies on using a proxy for 

the maximum circular v elocity of a halo dKlvpin. Truiillo-Gomez. & Prim 
l2010t|Prada et al.ll2.01 lb . The circular velocity is given by 

v 2 = GM(<r)/r (A2) 

For each halo in the sample we determine the maximum value 
of vl WK = max[GM/r] indirectly, by using the halo's mass pro- 
file to determine the RHS of Eq. IA21 and then divide by v^qo = 



16 



Dark Matter Halo Profiles of Massive Clusters 



i — i — i i ■ i ■ ■ i 1 — i — i i 1 1 1 ii 1 — i — i i 1 1 1 1 



T 



6 - 



full, A= 200 

• GS * HACC ■ G 

• G, center offset 



o 

O pL 
CM O 



4 - 



I 



U I I I I I I I I 



10 



13 



10 



14 



10 



15 



M 20 o h 1 M © 



6 - 



o 

O £- I 

CM O — 

O 



4 - 



~ i — i i 1 1 1 iij 1 — i i 1 1 1 iij 1 — i i ■ ■ ■ ■ ■ | 

\ full, A= 200 



(0-l-l)R v 




(O-l.OjRj 
--(0.1-2.0)R 
— (0.1-0. 5)R 

— (virialization) 

_i i i r 



vir 
vir 



_i i 



10 



13 



10 



14 



10 



15 



M 200 h 1 M © 



FlG. Al 1. — Tests to identify possible sources of systematic error. The left panel shows the three runs used in this paper and how the force resolution affects 
concentration measurements. There is excellent agreement between the GS and HACC runs, even though they were run using two completely different codes with 
different settings. The lower-resolution G run (resolution= 35/i~'kpc) is systematically lower by about 5% compared to the HACC ran (resolution= 7/r'kpc). We 
also study the possible effect of a misestimate of the halo center for the G run which can lead to a further reduction of concentration values. The right panel studies 
two other systematics issues, (i) the profile range used for fitting, and (ii) effect on concentration measurements when halos with large radial infall velocity are 
removed from the sample (solid curve). See the text for further discussion. 
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GMuoo/Rwo- Assuming an NFW form, one can relate vj^/vfoo 
to the concentration, C200, 

vjwc _ 0-2162C2QQ 
v| 00 m(c2oo) 

and solve for C2oo- The middle panel of Figure IA12I shows 
the ratio v max j V200 as a function of M200 for three redshifts for 
the relaxed cluster sample. The right panel compares the con- 
centrations obtained from the maximum velocity method and 
profile fitting. Again, the methods agree quite well within the 
error estimates, the velocity method being noisier. The mid- 
dle panel of Figure IA12I shows the ratio v max /v2oo as a func- 
tion of mass. Note that we cross over very smoothly from 
the GS run to the HACC run at M 2 oo = lO 14 /2 -1 M , more ev- 
idence for an excellent match between the results from these 
two simulations. One problem with this method is that because 
fmax ~ 2.2r s (where r max is the radius where v reaches v max ), 
at low concentrations, r maK becomes very close to r&, and is 
therefore very sensitive to any noise in the data, which will (i) 
result in biasing the concentration to a higher value (as seen in 
the flattening of the data at c ~ 3 in Fig. IA12I right panel), and 
(ii) make the result very sensitive to the shape of the measured 
profile at r ~ r&, increasing the possibility of system a tic er- 



rors . Finally, bo th Klypin. Truiillo-Gomez. & Primack I d2010t) 
and iPrada et al.l d201 lb have found that the maximum velocity 
method (along with their halo selection) leads to an upturn in 
the c-M relation at the high mass end at higher redshifts. We 
are unable to confirm this effect in our measurements, where 
we can investigate it (at redshifts, z= 1,2). 

We provide more information regarding the distribution of 
the concentrations, coqq, for the full halo sample at z = by con- 
sidering three representative mass bins, M200 = 5 x lO^/r'M©, 
1.5 x \Q w hT l M n and 8 x IO'VMq, as shown in Figure |3] 
Previous studies jJingl2000trShaw et al. 120061: iNeto et al. 120071: 
lDuffvetal.1l2008IK have fitted the concentration distribution to 
a log-normal distribution. However, this distrib ution is also 
very w ell described by a Gaussian as noted by iLukic et al. I 
( 120091) and lReed. Koushiappas. & Gap I d201 lb . Figure|3]shows 
that a Gaussian distribution provides a very good fit to our data, 
with relatively negligible non-Gaussian tails. 

We conclude that while statistical errors on the concentration- 



mass relation may have achieved ~ 5% accuracy in recent sim- 
ulation studies, systematic uncertainties of the order of 10% are 
apparently difficult to avoid. 
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